ON THE DEVELOPMENT OF BIOPHYSICAL MODELS FOR SPACE 
RADIATION RISK ASSESSMENT 

F. A. Cucinotta 1 and J. F. Dicello 2 


‘NASA, Johnson Space Center, Houston TX 77058-0001, USA 
2 The Oncology Center, Johns Hopkins University, Baltimore MD 21287-5001, USA 


ABSTRACT 

Experimental techniques in molecular biology are being applied to study biological risks from space 
radiation. The use of molecular assays presents a challenge to biophysical models which in the past have 
relied on descriptions of energy deposition and phenomenological treatments of repair. We describe a 
biochemical kinetics model of cell cycle control and DNA damage response proteins in order to model 
cellular responses to radiation exposures. Using models of cyclin-cdk, pRB, E2F’s, p53, and G1 inhibitors 
we show that simulations of cell cycle populations and G1 arrest can be described by our biochemical 
approach. We consider radiation damaged DNA as a substrate for signal transduction processes and 
consider a dose and dose-rate reduction effectiveness factor (DDREF) for protein expression. 

INTRODUCTION 

Theoretical models of initial damage to DNA have now reached a mature stage of development with 
Monte-Carlo track simulation codes that account for details of DNA structure, the hydration shell that 
interacts with DNA and most ionization and radical diffusion processes (Nikjoo et al., 1997; and Holley 
and Chatteijee, 1996). For describing cellular endpoints such as cell death, chromosome aberrations, or 
cell growth and arrest, models of the initial induction of lesions and phenomenological descriptions of 
lesion repair have been used to correlate experimental data for dose responses and the time development 
of cellular effects. However, everyone knows that the deleterious outcomes that follow radiation exposure 
are mediated and modulated by gene expression. Mathematical modeling of gene expression is a well- 
known approach in molecular biology, however has received little attention in the study of radiation 
effects Biochemical kinetics approach to modeling of mRNA and protein expression can be performed by 
computer simulation, including the solutions of non-linear differential equations that anse in a 
biochemical approach (Hargrove, 1995, Lauffenburger and Linderman, 1993). Important developments 
have been in studies of receptor binding and signaling (Lauffenburger and Linderman, 1993) and in 
studies of the behaviors of oscillatory enzymes (Goldbeter and Caplan, 1976). The role of diffusive 
processes in protein transport has also been considered (Goldbeter and Caplan, 1976; Lauffenburger and 
Linderman, 1993). The role of stochastic processes when low concentrations of molecular species are 
involved (Goss and Peccoud, 1998) may be important for signal transduction processes following low 
dose or low dose-rate exposures of cells. The mathematical description of covalent modification of 
enzymes (Goldbeter, 1981; Stadtman and Chock, 1977) predicted the importance of kinase activity in 
molecular controls related to cell cycle control and cancer including the large amplification of signals that 
may occur (Hunter, 1995). In this paper, we develop a mathematical description of several cellular 
processes using a biochemical model of the kinetics of cell cycle control and G1 checkpoint proteins. 


MATHEMATICAL MODEL OF CELL CYCLE CONTROL 


The pathways that control the cell cycle (Strauss et al., 1995) can mathematically be defined to 
successfully predict additional molecular interactions that modify the cell cycle. We have developed a 
model that incorporates sequential activation of cyclin-cdk complexes and their inhibitors. Identification 
of cyclin B and its associated cyclin dependent kinase cdkl (cdc2) led to several mathematical 
descriptions of the molecular interactions that control the G2/M transition (Tyson, 1991; Novak and 
Tyson, 1993). These authors differed in their treatments of the number of intermediate kinetic steps 
related to dimer formation between cyclins and cdk’s, the description of cyclin degradation, and control of 
cdk phophorylation sites. Two models of G1 control have been discussed (Obeyesekere et al., 1995, 
Hatzimanikatis et al., 1995) which have considered the pRb protein, however they did not consider other 
cyclins and cdks or their inhibitors. We describe the relationship of cdk4, cyclin D, pRb, E2F, downstream 
kinases, and cdk inhibitors such as pl6 and p21. We use the hypothesis that release of E2F/DP hetero- 
dimers from pRb provides a mechanism for the induction of downstream cyclins. We describe using non- 
linear differential equations, the periodic behavior and sequential activation of cyclm-cdks via coupled 
limit cycles (un-damped oscillatory solutions) and their transient behaviors. Our kinetic model uses 
compartments that describe the activation of each cyclin-cdk complex. Compartments are coupled by the 
release of the transcription factors E2F’s from pRb sequester upon phosphorylation of pRb by an upstream 
cyclin-cdk complex. Expression of cyclin D correlates directly with non-mutated pRb (Sherr and Roberts, 
1995) and we assume the cyclin D synthesis rate is proportional to the hypophosphorylated form of pRb in 
cells with normal pRb. Our model predicts that cyclin E is the primary regulator of cyclin A. Although 
this has not yet been observed directly, this conclusion is consistent with present experimental results. To 
complete the cell cycle and simulate its periodic behavior we require the dephosphorylation of pRb and 
the reformation of complexes of pRb with E2F s (Welch and Wang, 1995). 

The model kinetic equations for cyclin synthesis and degradation, cyclin-cdk interactions, transcription 
factor control, and cdk inhibitor induction can be written in a compact form as given in Table 1. Here we 
are using an indexing subscript that links proteins of a similar family or of a similar regulatory nature. In 
Table 1, k,j is rate constant for binding of molecules i and j, v s and v d (v P d) are synthesis and degradation 
rates, respectively. The k p (r p ) and k Q (r°) are rates of phosphorylation or de-phosphorylation, respectively. 
The k F is rate of re-formation of a pRb/E2F complex denoted [pRb-E2F]. Concentrations of molecules 
with ‘P’ superscript are activated form and with superscript ‘0’ inactive form. pRb is assumed to be 
dephosphorylated and to re-form dimers with E2F’s in G2/M as controlled by cyclin B-cdkl complex. The 
concentration of cyclin and cdk’s are denoted [Cycj] and [cdkj], respectively with the index i 
corresponding to cyclins D, E, A, etc. and the index j corresponding to cdk4, cdk2, etc. Transcription 
factors are denoted [TF,] and the transcription factors released by pRb are denoted [E2F,]. pRb, E2F, and 
cdk levels are assumed to be constant throughout the cell cycle (Welch and Wang et al. 1995). We do not 
consider individual members of the pRb family in the present model. In Table 1 the equations that control 
the time rate of change of phosphorylated pRb, denoted [pRb p ], and the free form of E2F are shown with 
the other forms of these molecules are solved by conservation. The formalism incorporates multiple 
inhibitors denoted as Inka for cdk4 and cdk6 inhibitors such as pl6 and pl8 (Lukas et al., 1995). 
Inhibitors bind specific cdk’s or cyclin/cdk complexes. Synthesis and degradation rates reflect signals 
from extra-cellular milieu. We have made a numerical study on cell progression by finding limit cycle 
solutions to the model. Values of rate parameters leading to periodic oscillations of protein concentrations 
are found by local stability analysis of the steady-state solutions of the model equations. Since we consider 
the sequential activation of 4 cyclin-cdk complexes, cell phase arrest corresponding to low or high cyclin- 
cdk complex activity at several points in the cell cycle are possible. Limit cycles found for the cyclin D- 
cdk4 complex determine the limit cycles of downstream kinases which are dependent on the release of 
E2F’s to switch from low activity modes. 


Table 1. Molecular Kinetics Equations for Cell Cycle Progression and Regulation. 


d[Cycj ]/dt = v Si [TFj] - v° Dl [Cyc,] - Zj k.j [Cyci][cdkj] 
d[Cyc P i ]/dt = Zj v Dlj [M P ij] - v p D , [Cyc p ,] 

d[cdkj]/dt = - Zi ky [Cyci][cdkj] + ZjVoij [M P j]+ ZjV°Di[M jj] - Z a k ja [cdkj][Inka] + 
£a V^a [Inka-cdkj] 

d[M°jj]/dt = kjj [Cycj][cdkj] - (k P jj+v°Di)[M 0 ,j] + k Q ij[M P j] - Zp kjjp[KIp][M y] + 

£pVDp[KIp-M°ij] 

d[M P ij]/dt = k P j[M°ij] - (v Di j + k Q y )[ M p y] 

d[pRb p ]/dt = Zj r P ij[M P ij][pRb-E2Fi+i]/(K+[pRb-E2Fi+i] - r Q G 2 /M [M p G2 /M][pRb P ] 
d[E2F i+ i]/dt = Zj r P ij[M P j][pRb-E2Fi+i] ]/(K+[pRb-E2F i+ i]) - k F i+i[E2F i+1 ][pRb 0 ] 
d[Inka]/dt = vsa [TF a ] - voa [Inka] - Sj kj a [cdkj][Inka] 

drinkq-cdkj]/dt = k ia [cdkjlflnkq ] - vpq pnkq-cdkj] 


In Figure la, we show the percentage of protein expression as a function of time for cychns E and A, and 
these results’ are compared with the data of Ohtusbo et al. (1995). Although they are not presented, the 
expression of free E2F’s and cyclin kinases will follow closely the peaks in cyclins D, E, and A. Figure lb 
shows the results of solutions for the expression of cyclin D, pl6 and the hyperphosphorylated form of 
pRb. These calculated results are compared with the data of Tam et al. (1994) for expression of pi 6. The 
differences between the predictions of our model and experimental measurements are largest after the S 
phase peak where the cell cultures become de-synchronized. Cyclin D is differentially expressed, while 
cyclins E, A, and B show similar behaviors in individual cell types. In Figure lc we show comparison of 
the model for cyclin E and cyclin A expression in HeLa cells (Pagano et al., 1992) which do not express 
cyclin D. The good agreement of calculations with the above experiments demonstrates that our model 
can describe the sequential activation of 4 distinct cyclins by using the controlled release of transcription 
factors from tumor suppressor proteins such as members of the pRb family. 

The growth kinetics of cell populations can be described by age-maturation diffusion equations mass- 
action rate equations, or discrete-time matrix methods. Here, if n,(t) denotes the number of cells in phase i, 
and the duration of the phase denoted as x„ the following set of differential equations describe the time- 
evolution of cell phase populations: 
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Figure la. Expression of cyclins E, A, and B in 
model and experiment for cyclin A and B (symbols). 


Figure 1b. Expression of cyclins D, p16, and pR 
in model and experiment for p16 (symbols). 
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Figure Id. Model prediction of cell phase 
kinetics for cells synchronized at M/G1 border. 


where the factor of 2 for the G1 population in equation (1) accounts of doubling at mitosis. The mass- 
action equations with constant rate coefficients will reflect exponential growth, however no fluctuations in 
individual phase populations will occur after about 2 mitosis which is inconsistent with normal cell 
growth. In order to relate molecular controls on cell cycle progression to cell population kinetics, we 
replace the constant coefficients in equation (1) with coefficients that are proportional to the activation of 
a cyclin-cdk complex for the Gl/S and G2/M transitions 

P p 

l/x s OC M E/cdk2, ; l/tG2 « M B/cdkl 


( 2 ) 










with results shown in Figure Id. This approach allows molecular switches described by our model to 
control progression of cell phase populations. This provides the mechanism for our model to directly 
relate protein fluctuations to cell phase population kinetics. Future applications will consider the loss of 
control of the mid-Gl restriction point and acceleration of G1 duration following oncogene activation or 
loss of tumor suppressor genes. 

P53 EXPRESSION AND SIGNAL TRANSDUCTION 

Early events in acute tissue responses are driven by cell killing and delays in cell division in highly 
proliferating cells. Cell killing can be separated into necrosis and apoptosis. These are distinguished by 
observations of swollen cell size, inflammation, and shrunken chromatin for necrosis and DNA laddering, 
fragmented nuclei, and lack of inflammation for apoptosis. It is well known that in most cell types the 
protein p53 controls the apoptotic response to radiation. Apoptotic cell death is prominent in many tissues 
important for acute health effects including lymphocytes, intestinal crypt cells, and bone marrow or colony 
forming units. Because of the high selection for p53 functional inactivation in human cancers, many cell 
culture models are deficient in apoptotic cell death. Dose-rate effects in apoptosis has been studied in a 
few cases and indicate less sparring with low LET radiation than what is normally observed for cell 
killing. Also, there are noted differences in molecular interactions of p53 and other molecules in important 
signal transduction pathways between human p53 and mouse p53 including the raf and MAPK proteins 
(Steegenga, et al 1996). The importance of such differences in radiation responses has not been 
investigated. Signal transduction following DNA damage by radiation is mediated by p53 coupling to the 
site of DNA breaks or DNA repair complex’s (Reed et al, 1995). Several outcomes are possible during 
signal transduction (Zhan et al., 1994; Canman et al, 1995) including growth arrest, programmed cell 
death (apoptosis) and also interference or cooperation with other signaling pathways as shown 
schematically in Fig. 2. It is believed that several of the p53 induced genes (PIGS) control the regulation 
of oxidative damage, ultimately leading to apoptosis (Polyak et al., 1997). The up-regulation of latent p53 
tetramers by post-transcriptional regulation is through the coupling of p53 tetramers to the DSB repair 
complex [C rep ] and converts it to the active form (Hupp and Lane, 1995). The activated form functions as 
a transcription factor. The PIGS are then regulated by /p53 le t] ■ 

We now discuss a mathematical model of the kinetics of p53 functions following DNA damage. Many of 
the details of these functions have not been established experimentally, and we make several simplifying 
assumptions in our model. First, because the tetramerization of p53 has been observed to be rapid relative 
to other cellular processes we consider only p53 tetramers. Second, following the ideas of Hupp and Lane 
(1994) we model the transactivation of p53 as an allosteric transition mediated by covalent modification 
of its carboxyl terminus and do not consider other phosphorylation sites in the present model. In many cell 
lineage’s, the presence of ssDNA or a DNA repair complex, leads to the activation of the trans-activating 
form of p53. We assume the presence of a DSB repair complex is the signal for activation of p53 in our 
present model. The latent form of p53 is known to undergo rapid turnover with a half-life of about 25 
minutes and the active form much longer-lived. The regulation of the transcription of p53 is modulated by 
the oncogene mdm2, which is transcribed following p53 activation. We write for the time-rate of change 
of the latent form of p53, 

= he,- v l d[ P 53 l ] - kA[signal][p53 L ] + kL[p53 A ] (3) 

dt 


and for the time-rate of change of the active form 



Fig. 2. Schematic of model reaction pathway of DNA damage leading to p53 expression, apoptosis or 
growth-arrest. 


d tP 5 ? A = kA[signal][p53 L ] -v A D [p53 A ]~ k L [p53 L ]~ kM[p53 A ][mdm2] (4) 

dt 

where the v D are degradation constants and the k are rate constants for (de)activation by conformation 
shifts of p53 or binding to mdm2 (k M ). Similar equations can be written for the transcription of mdm2 and 
binding of the mdm2 protein to activated p53. The binding of a single mdm2 protein to a p53 tetramer is 
believed not to block completely the tetramers function. We thus consider active forms of p53 that are 
complexed with one or more mdm2. We assume that p53 continues to function as a transcription factor, 
although at a reduced efficiency, until all 4 sites are bound by mdm2. We ignore any cooperative effects 
(positive or negative). Rate constants for mdm2 binding are then unchanged as each binding site becomes 
filled. The p53 transcription of other genes is reduced at the geometric rate as these sites are filled. We by- 
pass the intermediate step of the mdm2 RNA transcript that is easily included and should have only a 
small effect on modeling kinetics of the protein. 

The induction of DNA DSB’s is linear with dose for all dose-regimes of interest. Track structure studies 
indicate that the initial number of complex damages as defined by local damage clusters near the position 
of a DSB increase with LET or ionization density (Nikjoo et al., 1997, Holley and Chatteijee, 1996). 
Repair pathways are assumed to be dependent of the damage type which accounts for the slower kinetics 
observed in DSB with high LET radiation. The repair of DSB s of type j is described by 


d[DSB>] 


a 


dDose 


kl[DSBj][Erep ] 


( 5 ) 


dt 


dt 



Where a is the rate of formation of DSB’s (about 40 DSB’s per Gy of X-rays) and the time-rate of change 
of the repair-enzyme is described by 


d ^ ErSP l = -kl[DSB][Erep]+k2[Crep] (6) 

dt 

and of the repair-complex, 

d t CreP l = k\[DSB][Erep] - kj[Crep] ( 7 ) 

dt 

with the initial condition E 0 = [E rep J+[C rep ] . The reaction scheme described by equations (5)-(7) 
describes a dose dependent reaction velocity with saturation occurring at high doses (zero-order kinetics) 
and first-order kinetic at low doses. Studies of recombination repair of DSB’s indicate that there are 
several reaction pathways available including single strand annealing (SSA), and homologous and non- 
homologous recombination (Roth et ai, 1985; Fishman-Lobel et ai, 1992). We are studying the effects of 
repair pathway competition in our mathematical model. We also are exploring the role of base-sequence 
near the site of DBS’s which may effect processing such as SSA leading to deletions (Fishman-Lobel et 
al., 1992). We couple the [C rep ] to p53 as the ‘ signal ’ in equations (3) and (4). 

The kinase inhibitor p21 is normally found in quaternary complexes with the cyclins, cdk’s, and 
proliferating cell nuclear antigen (PCNA) (Namba et al. , 1995; Li et ai, 1994). Increases in p21 
concentration cause G1 arrest (Dulic et ai, 1994) through inhibiting G1 cyclin associated kinase activity. 
In our model we examine p53 dependent p21 upregulation after DNA damage as described by 

d[mRNA] = n + rp ^ p5 3 A ]- rD [ m RNA] ( 8 ) 

dt 

where r T and r D are the basal rates of transcription and degradation, respectively, and r p53 is the rate 
constant for coupling of the transcription factor p53 to the p21 promoter. The time rate of change of p21 
is given by 


V s[mRNA}-VD[p2l]-k,[p2l][M° D ]-ki[p2\][M 0 E ] (9) 

dt 

where v 5 and v D are the synthesis and degradation rates of p21, respectively. The last two terms in Eq. (9) 
are the coupling of p21 to the G1 cyclin-cdk complexes. For the background levels of p21 we use its 
known half-life of 30 minutes and assume that there is normally about one molecule of p21 in each of the 
G1 cyclin complexes. Fits of the model to the data of Bae et ai (1995) for the time course of p21 and p53 
expression in exposures of normal lymphoblast cells are shown in Figure 3. The comparison of the model 
to experiments at high doses determines model rate constants and allows for an alternate assessment of 
low dose-rate effects. In Figure 4 we show predictions of the model for the induction of DSB s, p53 and 
p2l expression at low dose-rate. These results indicate that although much sparring is seen in DNA 
damage repair, the persistence of p53 expression allows for a build-up of PIGS that would lead to 
apoptosis, growth arrest, and ultimately to harmful tissue effects. Modifications of the present model will 
be required to describe the modulation of DNA damage signaling proteins for time courses greater than 10 
hrs after DNA damage, including the role of phosphorylation or other conformation changes of p53 
throughout the cell cycle that are known to affect its DNA binding and transcription activity. 



Figure 3. Model calculations of p53 and p21 ex- Figure 4. Predictions of model for number of 

pression versus time after 6.3 Gy acute expo- DSB s remaining and relative expression of p53 

sure as fit to data of Bae et al, 1995. and p21 for 0.1 Gy/hr gamma irradiation. 

Genomic instability after radiation exposure would likely modify p53 levels because of the observed high 
sensitivity of p53 to DNA damage (Nelson and Kastan, 1994). These low dose-rate results predict that 
there is no threshold for p53 induction in agreement with observations of Lane (1998) for low dose-rate 
acute exposures. In Figures 5 we show comparisons of the model for acute exposures of 6.3 Gy in early 
Gl and at G1 + 6hrs. Figure 5a display a large G1 arrest in agreement with experiment (Dulic et al , 
1994). The results of the model are in general agreement with the Western blot analysis of Dulic et al. 
(1994). The results of Figure 5b show only minimal Gl arrest for the same exposure level. Here, the non- 
linear expression of kinase activity prevents the inhibitor in causing arrest. These results suggest that for 
asynchronous populations, Gl arrest is dependent on timing of signaling pathways and non-linear 
accumulation of cyclin associated kinases. 



Figure 5a. Calculations of relative expression 
of several proteins after acute exposure of 
6.3 Gy at the M/Gl border. 


Figure 5b. Same as Fig. 5a for exposure at 
6 hours past M/G 1 border. 




CONCLUSIONS 


Modification of intracellular signaling patterns and their effect on cell cycle control is a potential 
determinant in radiation response and must be adequately understood for extrapolating risk models to low 
dose-rates and the complicated radiation fields in space. We have described a mathematical model of cell 
cycle control through regulation of the related proteins and the coupling of this model to DNA damage 
through the p53 signal transduction pathway. Knowledge of protein-protein interactions, most 
importantly the residues involved in kinase activity and the possibly large number of binding proteins, 
may place limitations on models. However, biochemical models serve as a useful approach to consider 
these interactions and to model dose and dose-rate dependent responses. In future work, the model 
discussed here will be extended to describe exposures to asynchronous cell populations. The role of p53 
regulation by DNA damage many hours after exposure, and the description of the G2 arrest will be 
determinants in this description. A similar mathematical description can be applied to study other early 
events in the apoptosis pathway including p53 induction of the apoptosis inhibitor molecules, BAX and 
BAX formation of hetero-dimers with members of the apoptosis promoting Bel family. 
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